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■ We study the phase diagram at finite temperature of a system of Fermi particles on the sites of 

the Bethe lattice with coordination number z and interacting through onsite U and nearest-neighbor 
V interactions. This is a physical realization of the extended Hubbard model in the atomic limit. 
By using the equations of motion method, we exactly solve the model. For an attractive intersite 
' potential, we find, at half filling, a phase transition towards a broken particle-hole symmetry state. 

The critical temperature, as a function of the relevant parameters, has a re-entrant behavior as 
already observed in the equivalent spin-1 Ising model on the Bethe lattice. 
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I. INTRODUCTION 



U: 

Recently, it has been shown pj that a system of q species of Fermi particles, localized on the sites of a Bravais 
lattice, is exactly solvable in any dimension by means of the equations of motion approach [2| . This Fermi system is 
isomorphic to a spin-(j/2 Ising model in an external magnetic field. As a consequence, spin systems can be studied 
within a new approach 1] . Exactly solvable means that it is always possible to find a complete set of eigenvalues and 
'-^J . eigenoperators of the Hamiltonian, which closes the hierarchy of the equations of motion. Thus, one can get exact 
expressions for the relevant Green's functions and correlation functions. One finds that these functions depend on a 
finite set of parameters to be self-consistently determined [2| . It has been already shown how it is possible to fix such 
1^, ' parameters by means of symmetry and algebra constraints in the case of one dimensional systems for q = 1, 2, 3 [3[, 

and in the case of the Bethe lattice with coordination number z and q = 1 [j]. 
t-H ■ In this paper, we shall address the case of q = 2 species of Fermi particles on the sites of the Bethe lattice 
with coordination number z and interacting through nearest-neighbors interaction V . By considering also an onsite 
interaction U, one has the extended Hubbard model in the atomic limit. We find that for an attractive nearest- 
neighbor interaction there is a transition from a phase where the particle-hole symmetry is preserved to a phase 
where this symmetry is broken. The relative variations of the two regions depend on the coordination number z. 
Furthermore, for C//|l^| > 2 the phase diagram exhibits a re-entrant behavior. 

The plan of the paper is as follows: in Sec. [H] we outline the equations of motion method for the extended Hubbard 
model in the atomic limit. In Sec. IIII1 we compute the Green's and correlation functions and find that they depend 
only on two parameters. As a result, in Sec. IIV1 we are able to get a set of self-consistent equations enabling us 
• • . to fix these parameters in terms of which all the local properties of the system can be expressed. In Sec. [V] we 
shall determine the transition temperature as a function of the onsite interaction at half filling and analyze the phase 
diagram for a wide range of values of the parameters T, U (in units of \V\) and for different values of the coordination 
number z. In Appendix A we briefly discuss the equivalence between the extended Hubbard model in the atomic 
limit and the spin-1 Ising model on the Bethe lattice @, Q. Finally, Sec. IVIl is devoted to our concluding remarks. 



II. THE HAMILTONIAN AND THE EQUATIONS OF MOTION 



The extended Hubbard Hamiltonian in the atomic limit is given by 

h = + U T, D ® + \ v E (i) 

i i ~ i^j 

where U and V represent the onsite and the nearest-neighbor intersite interaction, respectively, fj, is the chemical 
potential, n(i) = n^ii) + n|(i) and D(i) — n^(i)ni(i) are the density and double occupancy operators at site i, 
respectively. As usual, n a (i) — cJ(i)C(j(«) with a — {1\ 1} and c a (i) (cj.(«)) is the fermionic annihilation (creation) 
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operator of an electron of spin a at site i, satisfying canonical anti-commutation relations. In the following we shall 
use the spinor notation for all fermionic operators. For instance, c^(z) = (c|(i) ct (*))- The sums in Eq. (TT|) run over 
the sites of a Bravais lattice. Here, we shall consider as Bravais lattice the Bethe lattice, which is an infinite Cayley 
tree consisting of a central site - which we denote by (0) - with z nearest-neighbors forming the first shell. Each site of 
a shell is joined to z — 1 nearest-neighbors to form the second shell, and so on to infinity. Thus, on the Bethe lattice, 
the Hamiltonian ^ may be conveniently written as 



H = -//n(0) + UD(0) + 2^ H {p > , (2) 
P =i 

where represents the Hamiltonian of the p-th sub-tree rooted at the central site (0) and it can be written as 

2-1 

H<*> = -mip) + UD(p) + Vn(0)n(p) + ^ H<*' m \ (3) 

m— 1 

where p are the nearest-neighbors of the site (0). In turns, ij(P' m ) describes the m-th sub-tree rooted at the site (p), 
and so on to infinity. 

The density operator cannot be used to employ the standard methods based on the equations of motion since it 
does not depend on time. In order to use the Green's function formalism one may consider the Hubbard operators 
= [n(i) — l]c(i) and r](i) = n(i)c(i), obeying to the following equations of motion: 

. d 



t^£(*) = -M£(i) + *^(i)n°(i) 

a 

i—r](i) — (U — jj)r]{i) + zVrj(i)n a (i) 



(4) 



where n a (i) = 2p=i n (hP)/ z i and (hP) are the nearest neighbors of site i. The algebra satisfied by the operators 
and D allows one to establish an important recurrence relation obeyed by n a (i): 



2z 

E 

m— 1 



A<$[n a (i)Y 



(5) 



(k) 

where the coefficients A m are rational numbers which can be easily determined by the algebra and the structure of 
the Bethe lattice, and which satisfy the relations J2m=i ^-m^ = 1 and Am = 5 m ,k {k = 1, • • ■ , 2z) @. The recurrence 
relation ([5]) is of seminal importance because it limits the number of composite operators which can be generated 
by the dynamics of the original fermionic operators. In fact, by taking successive time derivatives of the Hubbard 
operators and T](i), one clearly sees that for fc > 2z, no additional composite operators are generated and the 
equations of motion close Thus, one may define a new composite field operator ip(i) Q: 



( 



T](i) 
rj(i)[n a (i)] 



\ 



\v(i)[n a (i)? z J 



By means of the recursion formula (O, the operators and ip^^i) satisfy the equations of motion: 

. d 



»-^«(<) = [^«(i),fl] 



(6) 



(7) 



where e^' and e^' are the energy matrices of rank (2z + 1) x (2z + 1): 
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(8) 



(22 + 1) 
22 



3 



fU-fj, zV 
U -p, 







V o 






|(2*+1) 





zV 
U - n 






zVA\^ L > zVA^ z+1) 



U-p zV 

U-fJL 

zVA^+V zVA<i z z +p U-p + zVA^ 





zV 



(9) 



^22-1 



1 22 



whose eigenvalues, £^ and Em , are given by: 



-/i + (to - 1)F 

-At + tr + (m - i)v, 



(10) 



where m = 1, (2z + 1). It is easy to convince oneself that the Hamiltonian ((2]) has now been formally solved since 
one has a closed set of eigenoperators and eigenvalues. Then, by using the formalism of Green's functions (GF), one 
can proceed to the calculation of observable quantities. 

The two field operators ip^\i) and are decoupled at the level of equations of motion, as one may clearly see 

in Eq. ([7]). However, as we shall see in the next Sections, they are coupled via a set of self-consistent equations allowing 
for the determination of some unknown parameters in terms of which observables quantities may be computed. 



III. RETARDED GREEN'S FUNCTION AND CORRELATION FUNCTIONS 



The retarded thermal Green's function is defined as: 

G^(t-t') = (R [^)(o,i)^ t (o,o]) = ^-0({v^ ) (o,t),v (s)t (o,t')}), 



(11) 



where the index s refers either to the Hubbard operator £ or r\ and (• • ■ ) denotes the quantum-statistical average 
over the grand canonical ensemble. By means of the field equations ©, one finds that the retarded GF satisfies the 
equation 



w — e 



(12) 



where G^(u) is the Fourier transform of GW(< - if) and 1^ = (|-0 (j5) (O, i), ^ (s)t (0, is the (2z + 1) x (2z + 1) 
normalization matrix. The solution of Eq. (|12[) is : 



1 a (s,m) 



\ w-E$ +i6' 



(13) 



where Em) are the eigenvalues of the energy matrices, as given in Eq. (|10[). The spectral density matrices cr^ n> can 
be computed by means of the formula (2j: 



(s,n) 



(s,n) _ ( a ) 
'ab "an 



2z+l 



c=l 



l cb 



(14) 



In Eq. (fT4)) . f2( s ) is the (2z + 1) x (2z + 1) matrix whose columns are the eigenvectors of the energy matrix eW, It is 
straightforward to show that f2^ = fi^ = f2, where the matrix fl is given by: 



n 



p,k 



k = i,p = i 



I z \2z+l-y 

Moreover, the matrix elements of the normalization matrices in Eq. (|14[) can be cast in a simple form as [|[ 

T-(f) _ K (n+m-2) x (n+m-2) 
I n,m — K A 

t(v) — \(n+m-2) 



(15) 



(16) 
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where the charge correlators and A' fe ^ are defined as 

K « = (K(o)] fc ) 

( 

AW = i(n(0)K(0)] fe ). 
Similarly, one finds that the correlation function (CF) 

+00 

C {s) {t-t') = (V> w (0,t)V (s)t (0,O> = tt^t / ^e^^'Cl'ly, (18) 

—00 

can also be expressed in terms of the same charge correlators (fT7|) . In fact, by means of the relation 

f3uj~ 



1 + tanh : ^ 



/m [G(u>)] , (19) 



where /3 = 1/ &bT, the CF can be immediately computed from Eq. (TiU|) . One obtains 

2z+l 



m— 1 

22+1 



(20) 



C W(i -*') = - £ e- iE $ (*"'') a< s >") T«, 



2 

771=1 



where Ti s) = 1 + tanh(/3£^V 2 )- As one can clearly see, the knowledge of the GF's and, consequently of the CF's, 
is not fully achieved. In fact, they depend on the unknown correlators and A^which are expectation values of 
operators not belonging to the basis ©■ In the next Section we shall show how these quantities can be self-consistently 
computed. 

IV. SELF-CONSISTENT EQUATIONS 

In order to compute the unknown correlators (|17[) . one may start by noticing that, by exploiting the recursion relation 
([5]), also and A^ fe -* obey to a recursion relation which limits their computation just to the first 2z correlators 

2z 

( fc ) - \" 



K 

m— 1 

2z 



E 

(21) 

A« = £ 4 fc) A (m) . 

Then, one may split the Hamiltonian |(3J| as the sum of two terms: 

H = H ° +H > ( 22 ) 
Hi = zVn(0)n a (0). y ' 

Since Hq and Hi commute, the quantum statistical average of a generic operator O can be expressed as: 

_ Tr{Oe-^} (Oe-^)o 
W- Tr{e-? H } ~ (e-^)o ' 1 ' 

where (• • • )o stands for the trace with respect to the reduced Hamiltonian Hq. One then considers the correlation 
functions 

Ci ( S = (s(0)s t (0)[n Q (0)] ,! - 1 ), S = £,?7, k = 1, 2z + 1 , (24) 
which, by means of Eq. (|2li|) , can be written as: 
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The Pauli principle leads to the following algebraic relations 

? f (i)n(i) = rj f (i)n(i) = 

from which one has £^(0) e~^ Hl = £^(0) and rf(0) e~^ Hl = rj' (0) e - z VnC "(°) , The Hamiltonian i?o describes a system 
where the original lattice has been reduced to the central site (0) and to z unconnected sublattices. Thus, in the 
ffo-representation, the correlation functions connecting sites belonging to disconnected graphs can be decoupled. As 
a result, Eqs. ([25]) can be rewritten as: 



r ( S ) (£(0)e t (0))o(K(0)] fc - 1 )o r M (r / (0)ryt(Q)) o ([n"(0)] fc - 1 e-^" a W) o 



In the Ti/o-representation, the Hubbard operators obey to simple equations of motion: [£(i),ifo] = — /x£(i) and 
[r)(i), Hq] = — (/i — J7) r;(i). Thus, it is easy to show that the equal time CF's can be expressed as: 



/3m 1 



where: 



B 1 = {n(0)) = 26 (1 + C 



B 2 = <D(0)>o = 



1 + 2e^ + e /3(2/i-t/) 
1 + 2 e ^ + e^"* 7 ) 



and one has used the identities 



(5) _ l-j i +jg 

(B 1 -2g 2 )(e-^"°W) 
1,1 2<e-^> 

It is not difficult to show that the averages in the above equations can be expressed as: 

(e- pHl ) = 1-B 1 +B 2 + {B 1 - 2B 2 )(1 + aX x + a 2 X 2 ) z + B 2 (l + dX 1 + d 2 X 2 ) z 
(e-^ Vn °'^} = (l + aX 1 + a 2 X 2 y, 



(28) 



(29) 



Upon inserting Eqs. (|28|) into Eqs. (f27|) and by taking k = 1, one finds: 



(31) 



(32) 



where K — e ° v , a = (if — 1), d = (_ftT 2 — 1). Xi and X 2 are two parameters defined as: 

X x = (n Q ) , X 2 = (D a > . (33) 

Xi and X 2 are parameters of seminal importance since all correlators and fundamental properties of the system under 
study can be expressed in terms of them. Relevant physical quantities, such as the mean value of the particle density 
and doubly occupancy, and the charge correlators k 1 -^ and can be easily computed Q: 

(Xi - 2X 2 )(1 + aX x + a 2 X 2 ) + 2X 2 (1 + dX x + d 2 X 2 ) 



D = 



(1 - Xi + X 2 ) + (Xi - 2X 2 )(1 + aXi + a 2 X 2 ) + X 2 (l + dX t + d 2 X 2 ) ^ , . 

x 2 (i + dx 1 + d 2 x 2 ) ! ' ' 4 ' 

(1 - X x + X 2 ) + (Xi - 2X 2 )(1 + aXi + a 2 X 2 ) + X 2 (l + dXj + d 2 X 2 ) ' 
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and 

2z 



K (k) = 1 -^wry (K(°)] fc )o + E ( B i/™ z) + s 2ff W)(K(o)r+ fe ) 

A (fe) = 2{e - PHl)o { B i(K(°)] & )o + E + 2 ^)/i z) + 2B 2 ^)](K(0)] w+fe ) }> , 



m — 1 



/■(z) (z) 

where fm and c/m are some easily computable coefficients depending on the external parameters T and V and where 
the expectation value ([n a (0)] k )o can also be expressed in terms of X\ and X 2 Q. 

As a result, the solution of the model has been reduced to the determination of just two parameters. X\ and X2 can 
be determined by imposing translational invariance. In particular, the requests (n(0)) = (n a (0)) and (D(0)) = (D a (0)) 
lead to a set of two self-consistent equations Q : 

X 1 = 2e (3 »(l - Xi - dX 2 )(l + aXi + a 2 X 2 ) z - 1 + e J 3{2 ^- U) [2 + (d - 1)X 1 - 2dX 2 ] (1 + dX 1 + d 2 X 2 ) z ~ 1 

X 2 = e^- u ^ [1 + dX x - {2d + l)X 2 ](l + dXi + d 2 X 2 ) z - 1 - 2e^K 2 X 2 (l + aX x + c^Xtf- 1 . ' 

A full investigation of these equations will be given elsewhere. In this communication, we shall restrict our analysis 
to the possibility of a breakdown of the particle-hole symmetry. 

V. BREAKDOWN OF THE PARTICLE-HOLE SYMMETRY 

In this Section we shall study the possibility of a spontaneous breaking of the particle-hole symmetry. Let us 
consider the following value of the chemical potential: /i = U/2 + zV. For this value of /i, Eqs. ([55)) become 

X 1 = 2GK- z (l-X 1 -dX 2 )(l+aX 1 +a 2 X 2 ) z - 1 +K- 2z [2+(d-l)X 1 -2dX 2 }(l + dX 1 +d 2 X 2 ) z - 1 (37) 

X 2 = K- 2z [l + dX 1 -(2d+l)X 2 ]{l + dX 1 +d 2 X 2 ) z - 1 -2GK- z+2 X 2 (l + aX 1 +a 2 X 2 ) z -\ (38) 
and Eqs. (f3~4"|) become: 

2GK Z (1 + aX t + a 2 X 2 ) z + 2(1 + dX 1 + d 2 X 2 ) z 



D = 



K 2z + 2GK Z (\ + aXi + a 2 X 2 ) z + (1 + dX x + d 2 X 2 ) 

(1 + dx x + d 2 x 2 y 



(39) 



K 2z + 2GK Z (\ + aX x + a 2 X 2 ) z + (1 + dX x + d 2 X 2 ) z ' 
where G — e 131 1 1 2 . It is easy to show that for 

Xx = 1 - dX 2 (40) 
Eq. (|3"Tj) is always satisfied. Then, by substituting Eq. (|4"0"]) into Eq. (f3"9")l one obtains 

n = 1, D = — i — . — — . (41) 

2 + 2G(l-a 2 X 2 ) z V ' 

That is, solution (|40p is in agreement with the particle-hole symmetry. Upon inserting Eq. (|4"0")) into Eq. (f3"5)) one 
obtains: 

X 2 (l + K 2 ) + 2Gif A 2 (l - a 2 X 2 ) z - x -1 = 0. (42) 

Thus, for /i = U/2 + zV, Eqs. (f3T>)) admit a solution which satisfies the particle-hole symmetry and is described by 
the set of equations (f4"0|) and (|42[) . One may ask oneself if Eqs. ([57)) and (|55|) do admit other solutions, different from 
the one described by (|40|) and (|42|) and thus breaking the particle-hole symmetry. To this purpose, one may perturb 
the solution (|40p. by setting X% = 1 — e?X 2 + iu. With a little algebra it is easy to show that a solution with w ^ 
does exist if, and only if, the following equation 

2K Him+l ^ + R ^ _ 2 jj,_! + + _ l)*- 1 ^ - £T 2 (z - 2)] = (43) 
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is satisfied. This equation will fix the critical temperature T c below which the system admits other solutions, which do 
not satisfy the particle-hole symmetry. Numerical calculations show that Eq. (I43[) admits a solution only for negative 
values of V. Thus, in the following, we shall consider only the case of attractive intersite potential. It is interesting 
to consider Eq. (|4"3"| in two extremal limits, namely U / |V| — » — oo and T c — > 0. In the former limit, as expected, 
one finds that the critical temperature is the same as the one of a spinless fermionic system on the Bethe lattice, 
equivalent to the spin-1/2 Ising model [1, H, @j: 



2- 2 M ln(^) 



The limit T c — > provides the value of the ratio i//|V| at which one should expect a quantum phase transition: 



2e f3u/2 {z - 2) 2 - 1 - (z - \f- x K (z - 2) = => {U/2 - \V\) = k B T c \n 



(z-iy- 



2(z-2) 2 



(45) 



Thus, in the limit T c — > 0, one finds that = 2, independently of the value of z. The results obtained from Eq. 

(I43p are displayed in Fig. [1] One observes that, at fixed coordination number, by increasing U from large negative 
values (i.e., attractive onsite interaction) one finds a decrease of the critical temperature. 

An interesting feature of the model's phase diagram is that it shows regimes with a re-entrance: namely, fixing 
t//|V| > 2 and lowering the temperature, one switches from a particle-hole symmetry preserving phase to another 
one where this symmetry is broken at a certain critical temperature. Then, as evidenced in Fig. [U lowering further 
the temperature, one finds another critical temperature at which the symmetry is restored. How pronounced the 
re-entrance is, depends on z. 



VI. CONCLUDING REMARKS 



In this paper we have obtained the finite temperature phase diagram of a system of fermions with onsite and 
nearest-neighbor interactions localized on the sites of the Bethe lattice. The Hamiltonian describing such a system 
defines the so-called extended Hubbard model in the atomic limit. Upon using the equations of motion method, it 
is possible to exactly solve the model. For attractive nearest-neighbor interaction we find the critical temperature 
at which the system undergoes a transition to a phase where the particle-hole symmetry is broken. This critical 
temperature depends on the ratio U/\V\ and on z. It increases with increasing z and presents a re-entrant behavior 
for U/\V\ > 2. 



8 



APPENDIX A: SPIN-1 ISING MODEL ON THE BETHE LATTICE 



In this Appendix we shall analyze the correspondence between the extended Hubbard model on the Bethe lattice 
and the spin-1 Ising model defined on the same lattice. A transformation from a fermionic to a spin Hamiltonian can 
be performed by the use of the pseudospin variable S(i): 

S(i) = ni(i)+n l (i)-l = n(i)-l. (Al) 

S(i) can take four values, with S(i)=0 double degenerate: 



(A2) 



Under the transformation 



n r (i) = 


= 


- 




S(i) = 


-1 


n T (i) = 1 


n l(i) = 


= 




S(i) = 





n T (i) = 


n l(i) = 


: 1 




s® = 





n T (i) = 1 


n l(i) = 


: 1 


<*=> 


S(i) = 


1 


the Hamiltonian ([2]) can 


be cast in 


the form: 




H = AS 2 (0) 


-hS(0)- 


f Eq - 


z 

^E 


H (p) 





ffW = -JS(0)S(p) + A S 2 ( P ) ~hS{ P )+ H 



(A3) 



m— 1 



where E = -V + (-fi + V)N, N is the total number of sites, h = fx - zV - U/2, J = -V and A = U/2. Thus, 
the Hamiltonian (|A3|1 appears as the one of a spin-1 Ising model with nearest-neighbor exchange interaction J in the 
presence of a crystal field A and of an external magnetic field h. The difference here is that the Hamiltonian (|A3[) 
is pertinent to a four-level system because of the spin degeneracy. It is possible to get rid of the spin degeneracy 
by mapping the fermionic Hamiltonian on the standard spin-1 Ising one with S(i) = { — 1,0, 1} paying the price of 
making the crystal field A to be temperature dependent 0,0): A = U/2 + fc^Tlog 2. The double degeneracy of every 
S(i) — leads to a factor 2 for every singly occupied site in the partition function of the classical spin system. This 
gives rise to an overall factor 

JJ 2 i-s 2 M =2 E<[i-s 3 W] ( A4 ) 



One may rewrite the partition function of Hamiltonian (|A3|) as follows: 

Z= J2 exp{-(3H[S(i)}} = E exp{-(3H[S(i)}} (A5) 

{S(<)}={-1,0,0,1} {S(i)}={-l,0,0,l} 

where H is the Hamiltonian of the standard spin-1 Ising model on the Bethe lattice, but now with an effective 
temperature-dependent crystal field: 



H = AS 2 (0) — h S(0) + E + ^2H 



(p) 



(A6) 



z-l 



H {p) = -JS(0)S(p) + A S 2 (p) - h S{ P ) 



m— 1 



where Eq = Eq + fcsTln2 and A = A + ksT In 2. Having established the mapping between the two models, we find 
that our critical temperature exactly agrees with the one previously found in the literature @. 
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